use std::f32::consts::E;
use std::f64::consts::PI;

// 欧拉-马歇罗尼常数
const EULER_MASCHERONI: f64 = 0.57721566490153286060651209;

fn sigmoid(x: f64) -> f64 {
    return 1.0 / (1.0 + ((-x).exp()));
}

fn softplus(x: f64) -> f64 {
    return (1.0 + (x.exp())).ln();
}

fn relu(x: f64) -> f64 {
    return if x >= 0.0 { 1.0 } else { 0.0 };
}

fn factorial(n: f64) -> f64 {
    if n == 0.0 {
        return 1.0;
    }

    return n * factorial(n - 1.0);
}

fn gamma(x: f64) -> f64 {
    if x == 1.0 {
        return 1.0;
    } else if x == 0.5 {
        return PI.sqrt();
    } else if x >= 30.0 {
        return stirling(x);
    } else if x > 1.0 {
        return (x - 1.0) * gamma(x - 1.0);
    } else if x < 0.0 {
        return gamma(x + 1.0) / x;
    }

    //魏尔斯特拉斯无穷乘积
    let mut m = x * ((E as f64).powf(EULER_MASCHERONI * x));
    for n in 1..=100 {
        m *= (1.0 + (x / n as f64)) * ((E as f64).powf(-x / n as f64))
    }

    return 1.0 / m;
}

fn beta(m: f64, n: f64) -> f64 {
    gamma(m) * gamma(n) / gamma(m + n)
}

fn stirling(x: f64) -> f64 {
    (2.0 * PI * x).powf(0.5) * ((x / (E as f64)).powf(x))
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_gamma() {
        println!("{:?}", factorial(60.0));
        println!("{:?}", gamma(60.0));
    }
    #[test]
    fn test_beta() {
        println!("{:?}", beta(60.0, 30.0));
    }

    #[test]
    fn test_relu() {
        println!("{:?}", relu(10.0));
        println!("{:?}", relu(-10.0));
    }

    #[test]
    fn test_softplus() {
        println!("{:?}", softplus(10.0));
    }

    #[test]
    fn test_sigmoid() {
        println!("{:?}", sigmoid(10.0));
    }
}
